The role of conventional and unconventional adaptive routes in lowering of intraocular pressure: Theoretical model and simulation

In this article, we propose a theoretical model leveraging the analogy between fluid and electric variables to investigate the relation among aqueous humor (AH) circulation and drainage and intraocular pressure (IOP), the principal established risk factor of severe neuropathologies of the optic nerve such as glaucoma. IOP is the steady-state result of the balance among AH secretion (AHs), circulation (AHc), and drainage (AHd). AHs are modeled as a given volumetric flow rate electrically corresponding to an input current source. AHc is modeled by the series of two linear hydraulic conductances (HCs) representing the posterior and anterior chambers. AHd is modeled by the parallel of three HCs: a linear HC for the conventional adaptive route (ConvAR), a nonlinear HC for the hydraulic component of the unconventional adaptive route (UncAR), and a nonlinear HC for the drug-dependent component of the UncAR. The proposed model is implemented in a computational virtual laboratory to study the value attained by the IOP under physiological and pathological conditions. Simulation results (i) confirm the conjecture that the UncAR acts as a relief valve under pathological conditions, (ii) indicate that the drug-dependent AR is the major opponent to IOP increase in the case of elevated trabecular meshwork resistance, and (iii) support the use of the model as a quantitative tool to complement in vivo studies and help design and optimize medications for ocular diseases.


I. INTRODUCTION
In this article, we propose a theoretical model to investigate the processes of aqueous humor (AH) circulation and drainage in the human eye [1][2][3][4] under the assumption of a given volumetric flow rate (VFR) of AH secreted by the ciliary processes.
The main clinical motivation of this study is that AH establishes the steady-state value of the intraocular pressure (IOP) in the globe of the eye. 5Lowering of IOP is the only currently approved therapeutic approach for a severe eye disease such as glaucoma, so that, being able to keep IOP at a controlled, healthy level, typically around 15 mm Hg, may prevent retinal ganglion cells and optic nerve head from irreparable damage. 6,7A pressure level of 21 mm Hg is usually taken as a measure for screening or diagnosis of glaucoma, even if several population-based studies indicate a great interindividual variation in the susceptibility of the optic nerve to IOP-related damage. 8n abnormal value of IOP is the macroscopic result of multiple disorders localized at very different spatial scales, ranging from 1 nm ¼ 10 À9 m (level of cellular plasma membrane of the ciliary epithelium) to 1 mm ¼ 10 À3 m (level of the trabecular meshwork, TM).
To deal with such a wide spatial range in a computationally effective manner, in the present work, we develop a mathematical model of AH dynamics based on the average of the laws of mass and linear momentum balance of the aqueous humor fluid over the posterior chamber (PC), the anterior chamber, the trabecular meshwork, and the ciliary muscle.The result of the average is a zero-dimensional (0D, or compartment) representation of the anterior segment of the eye that includes • a compartment representation of AH circulation from the posterior chamber (PC) into the anterior chamber (AC); • a compartment representation of AH drainage through trabecular and uveoscleral (US) outflow pathways.
The advantage of averaging mass and momentum balance laws over the various regions of the anterior segment of the eye is that the resulting compartment representation can be conveniently described in mathematical terms by leveraging the analogy between fluid and electric variables illustrated in Ref. 9. AH circulation is the result of passive fluid diffusion driven by hydraulic pressure gradients across the series of two linear hydraulic conductances representing the posterior and anterior chamber. 10AH drainage is the result of adaptive mechanisms modeled by the parallel of two main outflow pathways, a "conventional" adaptive route (ConvAR) throughout the TM 11,12 and an "unconventional" adaptive route (UncAR) throughout the ciliary muscle. 13A significant novelty of the formulation proposed in this article compared to the existing literature is the mathematical treatment of the UncAR.Based on the studies of Johnson and Erickson, 14 Johnson et al., 15 and Costagliola et al., 13 the UncAR is modeled by the parallel of a pressure-dependent route and a drug-dependent route.The former is described by a nonlinear hydraulic conductance depending on the pressure difference between the AC and the suprachoroidal (Sch) space, 16 the latter by a current source nonlinearly controlled by the amount of topical drug administered to lower the IOP.
The model is numerically implemented in a computational virtual laboratory (CVL) for subsequent use in the study of problems of clinical interest.In particular, the focus is devoted to determine the value attained by IOP under physiological and pathological conditions.In the former case, model parameters are calibrated to reproduce the characteristic values of AH volumetric flow rate and IOP in a normal healthy eye. 3 In the latter case, pathological conditions are simulated by increasing the value of the TM resistance with respect to baseline conditions. 17,18Results from extensive simulations with the CVL: (1) confirm the conjecture of Costagliola et al. 13 that the UncAR acts as a relief valve under pathological conditions related to an abnormal increase of TM resistance; (2) indicate that the drug-dependent AR is the major opponent to IOP increase due to elevated TM resistance; (3) support the use of the model as a quantitative tool to complement in vivo studies and help design and optimize medications for ocular diseases.
An outline of the article is as follows.Section II gives an introduction to the physiological processes of AH dynamics in the human eye.Section III introduces the three-dimensional (3D) geometric and differential model to study AH circulation and drainage in the human eye.Section IV illustrates the procedure to reduce the 3D geometric and differential model to a compartment model based on a suitable average of the mass and momentum balance equation over each region of the anterior segment of the eye.In Sec.V, the reduced model is translated into a conceptual model by leveraging the analogy between fluid and electric variables, while in Sec.VI, the conceptual model is translated into a system of algebraic equations expressing Kirchhoff's law for the current at each node of the circuit and Ohm's law for the linear and nonlinear hydraulic conductances.Section VII describes the numerical values of model parameters that are used in the numerical simulations, and Sec.VIII describes the algorithm, based on a fixed-point iteration, to compute the intraocular pressure as a function of the other input variables and parameters.Section IX contains a detailed report of the computational tests that are performed to validate the proposed model and fixed-iteration algorithm in the study of two clinically relevant situations.Section X discusses the physiological consistency of the simulation results illustrated in Sec.IX and the clinical relevance of the adaptive routes of the unconventional outflow pathway on the lowering of IOP.Section XI summarizes the main findings of the article, indicating some limitations of the proposed model and suggestions to extend its use as a flexible tool to complement clinical practice and experimental activity.Section XII collects the figures illustrating the simulation results described in Sec.IX and discussed in Sec.X. Appendixes A and B illustrate the main results of the theoretical analysis of the solution map, providing sufficient conditions to be satisfied by model parameters to ensure the existence and uniqueness of the solution and the convergence of the fixed-point iteration.

II. PHYSIOLOGY OF AQUEOUS HUMOR DYNAMICS
Figure 1 illustrates the processes contributing to AH dynamics, whereas Fig. 2 illustrates the structure of the ciliary epithelium.Aqueous humor is secreted by the ciliary epithelium into the PC, which is the region behind the iris and in front of the ciliary body and the lens.Then, the AH fluid is not stagnant, rather, it undergoes the following physical processes that regulate its motion throughout the anterior segment of the eye: M.1: AH circulation from the PC, through the pupil, to the iridocorneal angle of the AC, shown in Fig. 3; The TM outflow pathway is referred to as the conventional adaptive route (ConvAR) and includes the TM, the Schlemm canal, the collector channels, and the episcleral venous system.The US outflow pathway is referred to as the unconventional adaptive route (UncAR) and includes the iris root and the interstitial spaces of the ciliary muscle through which AH flows into the suprachoroidal space (see Fig. 4).The motion of the AH fluid through the anterior segment of the eye is driven by • the hydraulic pressure difference between PC and AC, AC and episcleral vein, and AC and Sch space; • the temperature difference between the iris surface (37 C) and cornea surface (25 C).
The hydraulic contribution determines the magnitude of AH volumetric flow rate, whereas the thermal contribution controls the AH direction across the AC.Aqueous humor viscosity may be an important factor in AH flow of glaucomatous eyes, 19 as can also be seen from relations (6) and ( 7) which show that, for a given pressure drop, the aqueous fluid velocity increases as the fluid becomes more inviscid and decreases if the fluid gets more viscous, as expected from physical intuition.Fluid temperature is responsible for the fact that AH flow is almost tangential to the iris surface, as depicted in Fig. 1.This is numerically confirmed by the 3D finite element simulations by Wang et al. 20 where the mass and momentum balance equations for the fluid are coupled to the energy balance equation for the temperature in which the Boussinesq approximation is adopted to relate fluid density to temperature.

III. 3D MODEL FOR AQUEOUS HUMOR DYNAMICS
In this section, we formulate the three-dimensional (3D) geometric and differential model to study the processes of AH circulation and drainage in the human eye.

A. 3D geometric representation
Figure 5 shows a cross-sectional view of the geometrical representation of the anterior segment of the eye considered in this article.The PC and AC are two disks obtained by rotation around the pupillary axis of the green and magenta rectangles, respectively.The pupil is a cylinder obtained by rotation around the pupillary axis of the gold rectangle.The iris is a hollow disk obtained by rotation around the pupillary axis of the brown rectangles, which separates the PC disk from the AC disk.The PC disk is connected to the ciliary capillary pressure p cc by the ciliary epithelium (CE) compartment block.The red arrows indicate the volumetric flow rate Q in that is secreted by the ciliary epithelium and enters the PC.The cyan and yellow rectangles are subsets of the PC and AC, respectively, and represent the transition regions for fluid passage from the PC, through the pupil, into the AC, as indicated by the thick black arrows crossing the vertical dashed lines.The dark blue arrows indicate the AH volumetric flow rate Q out that is drained out of the AC.This flow rate is divided into two components.The first component of Q out goes into the ConvAR, which is connected to ground by the episcleral vein pressure p ev .The second component of Q out goes into the UncAR, which is connected to ground by the suprachoroidal space pressure p S .

B. 3D fluid model
Let us denote by X PC ; X PC;tr ; X IOP , v X AC;tr , and X AC the 3D subdomains that are obtained by rotation of the green rectangles; the cyan, gold, and yellow rectangles; and the magenta rectangles in Fig. 5 around the pupillary axis, and let The surface delimiting the boundary of X is denoted by R, and the outward unit normal vector on R is indicated by n.The boundary surface R is the union of the disjoint boundaries R in (interface between the PC and the CE which secretes the AH flowing into the PC), R len (interface between the PC and the lens), R int (internal surface separating the PC from the iris, the iris from the pupil and the iris from the AC), R cor (interface between the AC and the cornea), and R out (interface between the AC and the outflow pathway region, comprising the trabecular meshwork and the ciliary muscle).
Throughout the remainder of the article, we make the following assumptions: (a) constant fluid density (water), (b) no net production of fluid inside X, (c) stationary flow, (d) small Reynolds number, and (e) no gravitational effects.Under assumptions (a)-(e), the 3D model for the AH fluid is represented by the following differential system for fluid velocity v and fluid pressure p: where r s v ¼ ðrv þ ðrvÞ T Þ=2 is the symmetric gradient of v (the strain rate, units: s À1 ) and l fl is the dynamic viscosity of the fluid (units: Pa s).The system of partial differential equations (PDE) (1) constitutes the conservative form of the Stokes model for a viscous incompressible fluid in stationary motion in the 3D domain X under the action of a given volumetric flow rate Q and given external pressures p cc ; p ev , and p s .The boundary conditions that complement system (1) read on R in : Equation (2a) expresses the balance between the volumetric flow rate that flows into the PC and the volumetric flow rate that flows out of the PC.Equation (2b) is the functional relation between the trace p PC of the fluid pressure on R in on the side of the PC and the blood pressure p cc in the ciliary capillaries.Equation (2c) expresses the no-slip condition for the aqueous fluid at the wall boundaries.Equation (2d) is the functional relation between the trace p AC of the fluid pressure on R out on the side of the AC and the episcleral vein pressure p ev and the suprachoroidal space pressure p S .

IV. REDUCED MODEL FOR AH FLOW
The solution of the 3D continuum model ( 1) and ( 2) requires an intensive computational effort because of the coupling between the description of AH motion inside the domain X and the description of AH secretion, represented by the functional relation (2b), and AH drainage, represented by the functional relation (2e).As the focus of this work is on the study of the role of conventional and unconventional outflow pathways in the drainage of AH, an alternate approach based on model reduction is desirable, especially in the case where multiple simulation tests have to be performed on different input datasets.

A. Reduction of momentum balance
The first step to derive the reduced model for AH flow in the anterior segment of the eye is to solve the momentum balance equation (1b) under suitable simplifying assumptions.
To this purpose, we introduce the two-dimensional domain X 2D illustrated in Fig. 6 to represent any of the green or magenta rectangles.We denote by x, y the spatial coordinates in X 2D , and e x and e y the unit vectors of the x and y axes, respectively.Then, we make the following assumptions: A.1: The spatial distributions of fluid velocity and pressure in the 3D domain X do not depend on the rotation angle around the pupillary axis.A.2: The fluid velocity vector v has the profile of Poiseuille flow along the x coordinate, i.e., v y ðx; yÞ ¼ 0 8ðx; yÞ 2 X 2D ; (3a) where Ṽ is a constant and uðyÞ ¼ 3 2 ð1 À ðy=aÞ 2 Þ, in such a way that Ṽ coincides with the mean integral value of v x over ½Àa; a. FIG. 5. Cross-sectional view of the geometrical scheme of the anterior segment of the eye.The green, gold, brown, and magenta rectangles represent the posterior chamber, the pupil, the iris, and the anterior chamber, respectively.The other symbols and objects and their meaning are described in the main text.
FIG. 6.The domain X 2D represents any of the green or magenta rectangles in Fig. 5.The length is b and the height is 2a.

Using assumptions A.1 and A.2 into the constitutive law (1c) for the fluid stress tensor yields
Tðx;yÞ ¼ 2l fl @v x ðx;yÞ @x À pðx;yÞ l fl 2 @v y ðx;yÞ @x þ @v x ðx;yÞ @y !l fl 2 @v y ðx;yÞ @x þ @v x ðx;yÞ @y !2l fl @v y ðx;yÞ @y À pðx;yÞ Replacing the above expression into (1b) yields The second equation in (4) implies that p does not depend on y, while the first equation in (4) yields Since the left-hand side in ( 5) is a constant, we see that the fluid pressure varies linearly with x, and the velocity in the PC and AC is given by the following discrete Darcy's law: where Dp :¼ pð0Þ À pðbÞ and is the specific hydraulic conductance of the 2D domain X 2D (units: m s À1 Pa À1 ).

B. Reduction of mass balance
The second step to derive the reduced model for AH flow in the anterior segment of the eye is to solve the mass balance equation (1a) under suitable simplifying assumptions.
To this purpose, we denote by V a bounded open domain of R 3 with sufficiently smooth boundary @V and we denote by n @V the unit outward normal vector on @V.For any differentiable vector field Z : V !R 3 (units: Z), we denote by the net flux of Z across the volume boundary @V (units: Z m 2 ).In the case where Z is a velocity field (units: m s À1 ), the next flux U V ðZÞ has the physical meaning of volumetric flow rate (units: m 3 s À1 .Then, we introduce four fluid parallelepipeds, P PC ; P AC , P Conv , and P Unc , and one node, P IOP , illustrated in Fig. 7 to represent the compartment equivalent scheme of the anterior segment of the eye.
We let x PC ¼ X PC [ X PC;tr and x AC ¼ X AC [ X AC;tr , and we denote by X TM and X cm the three-dimensional regions of the anterior segment of the eye that are occupied by the TM and the ciliary muscle, respectively.The parallelepipeds and the node enjoy the following properties: P.1: U PUnc ðvÞ ¼ U Xcm ðvÞ: • P.2: the volume of X IOP is "lumped" into the node P IOP ; • P.3: the aqueous fluid flows only through the upper and lower surfaces of each parallelepiped.Property P.1 expresses the fact that each parallelepiped is equivalent to its ocular counterpart from the point of view of mass balance.Property P.2 expresses the fact that the volume of the pupil is neglected with respect to the volume of the other compartments.Property P.3 expresses the fact that the aqueous fluid is assumed to flow in the direction connecting the top and bottom surfaces, as represented in Figs. 6 and 7.
We make the following assumptions: • P.4: The fluid pressure in P PC is constant and equal to p PC (indicated by the red bullet in Fig. 7); • P.5: the fluid pressure of node P IOP is constant and equal to IOP (indicated by the black bullet in Fig. 7); • P.6: the fluid pressure in P AC is constant and equal to p AC (indicated by the dark blue bullet in Fig. 7); P.7: the functional relation (2b) is defined as where Q is a given input parameter representing the volumetric flow rate secreted by the CE P.8: the functional relation (2e) is defined as where Q Conv represents the AH volumetric flow rate drained out through the ConvAR and Q Unc represents the AH volumetric flow rate drained out through the UncAR.The mathematical expressions of Q Conv and Q Unc as a function of the biophysical properties of P Conv and P Unc will be provided in Sec.VI.We now write the mass balance equation (1a) over x PC , X IOP , and x AC .Using property P.3, P.2, and P.7, in conjuction with the discrete Darcy law (6), we obtain the reduced-order model of AH circulation and drainage The quantities L xPC and L xAC are the hydraulic conductances of the PC and AC, respectively (units: m 3 s À1 Pa À1 ), defined as R PC and R AC being the cross-sectional areas illustrated in Fig. 7.The cross-sectional areas cannot be explicitly characterized as a function of measurable biophysical parameters so that the hydraulic conductances must be determined through a calibration procedure, as described in Sec.VII.The quantities Q Conv and Q Unc are the volumetric flow rates of the conventional and unconventional adaptive routes.Their mathematical characterization will be provided in Sec.VI.

V. THE CONCEPTUAL MODEL
In this section, we translate the reduced model ( 11) into an equivalent electric circuit.The adoption of this approach in ocular pathophysiology is well known 2 and has a twofold advantage: (a) it allows us to write AH model equations in terms of Kirchhoff's current law, expressing mass conservation at each node of the circuit, and Ohm's law for any hydraulic conductance in the circuit, expressing the discrete Darcy law (6), and (b) it helps us explicitly characterize Q Conv and Q Unc by means of their electric analogs.
Figure 8 is a conceptual model of the various processes described in Sec.II and is the translation of the reduced-order model (11) based on the analogy between fluid variables and electric variables described in Sacco et al. 9 according to which electric potential (units: V) is the analog of fluid pressure (units: mm Hg), while electric current (units: A) is the analog of volumetric flow rate (units: m 3 s À1 ).
The input data of the model are the volumetric flow rate Q of AH, represented by the current source at the top of the scheme in Fig. 8, and the voltage sources p ev and p S , representing the blood pressure in the episcleral vein (ev) and Sch space, respectively.The current source is the equivalent electric parameter that is used in the present article to describe AH secretion from the ciliary processes of the eye.We refer to Raviola and Raviola 21 for a description of the cellular and subcellular structure of the ciliary epithelium (CE), Sala et al. 22 for a mathematical model and simulation of transmembrane fluid flow, and Sacco et al. 23 for a mathematical model and simulation of the AH secretion process.
The upper block of the scheme represents the process of AH circulation between the PC and the AC.The symbols L PC and L AC represent the hydraulic conductance of the PC and AC, respectively (units: m 3 s À1 mm Hg À1 ); the quantities p PC and p AC represent the pressure exerted by the fluid in the PC and AC, respectively (units: mm Hg), while the quantity IOP represents the intraocular pressure resulting from the balance among secretion, circulation, and drainage (units: mm Hg).
The lower block of the scheme represents the process of AH drainage.The quantities Q Conv and Q Unc are the volumetric flow rates of AH leaving the eye through the conventional and unconventional pathways, respectively (units: m 3 s À1 ).The quantity L Conv is the hydraulic conductance of the TM (units: m 3 s À1 mm Hg À1 ).The

VI. MODEL EQUATIONS
The application of Kirchhoff current law (KCL) at nodes 1, 2, 3, and 4 in the electric scheme of Fig. 8 leads to the following system of algebraic equations: The quantities Q The application of the discrete Darcy law ( 6) combined with the definition of the hydraulic conductances (11d) and (11e) leads to the following set of constitutive relations: Replacing relations ( 13) into (12) yields where m (units: g) is the strength of topical drug in a given dosage form (1 drop/day), defined as the drug concentration times the volume of one drop of the drug.Relations (14) constitute a system of nonlinear algebraic equations for the dependent variables x, y, z, and w.The system can be reduced into a single nonlinear algebraic equation for the sole variable x by considering y, z, and w as inputs and expressing them as a function of x.Replacing Eq. (14a) into (14b) yields where R AC ¼ L À1 AC is the hydraulic resistance of the AC (units: mm Hg s m À3 ).Using now (14c) and ( 15), we get where Conv is the hydraulic resistance of the TM (units: mm Hg s m À3 ) and Q Unc can be computed from (14d).
Relation ( 16) is a mathematical translation of the conjecture that the UncAR acts as a relief valve to aqueous humor water flow under pathological conditions. 13According to the second term on the righthand side of ( 16), an increase of TM resistance, R Conv , turns out into an increase of intraocular pressure, which is counterbalanced by the negative pressure drop ÀR Conv Q Unc associated with the UncAR.
To solve Eq. ( 16), we need to provide a functional relation for L hyd Unc (with respect to the system unknown w) and Q drug Unc (with respect to the input datum m).With this aim, we refer to Keener and Sneyd 24 and introduce the Hill function associated with the variable n !0, where p is the Hill coefficient, while b and K act are positive parameters.
From now on, in each occurrence, we set p ¼ 4. We see that f H is a sigmoidal monotonically increasing function of n such that f ð0; b; p; The rapidity at which f H reaches the asymptotic value b is determined by K act , see Fig. 9 in the case b ¼ 1.
For a generic quantity U, we denote by U b the value of U in baseline conditions, and we introduce the relationships In the remainder of the text, we assume that the AH volumetric flow rate Q, the suprachoroidal pressure p S , and the AC resistance R AC are set equal to their baseline values Q b ; p S;b , and R AC;b , respectively.This assumption implies that where Let us characterize the values of the parameters b max and K act for the pressure-and drug-dependent components of the UncAR.In the case of the hydraulic conductance L hyd Unc , we set In the case of the volumetric flow rate Q drug Unc , we evaluate b drug unc;max using relation (B2e).The parameter K act;drug is given by 0:5M tot;abs , where M tot;abs (units: lg) is the total mass of topical drug that is actually absorbed by the anterior segment of the eye, given by the following expression: where F and m onedrop denote the bioavailability of Latanoprost 0.005% and the mass contained in one drop (units: lg), respectively.The corresponding values of the two parameters are b drug unc;max ¼ 7:87 Â 10 À11 m 3 s À1 ; K act;drug ¼ 0:547 lg: The graphs of ZðxÞ and the corresponding hydraulic conductance L hyd Unc ðZðxÞÞ as a function of x (intraocular pressure) varying between 10 and 30 mm Hg are illustrated in Fig. 10.
Equations (20) are phenomenological relations describing a saturation trend as a function of an increasing parameter.In the case of the pressure-dependent component of the unconventional outflow pathway, the increasing parameter is the deviation from baseline conditions of the pressure difference between anterior chamber and suprachoroidal space, consistent with physiology evidence shown in Refs.13 and 16.In the case of the drug-dependent component of the unconventional outflow pathway, the increasing parameter is the amount of drug that is administered in a therapy to cure an ocular disease, consistent with the observation that arbitrary increasing drug concentration does not significantly improve the efficacy of the adaptive response of the system.Based on these considerations, the Hill function (17)  appears to be a good candidate represent the saturation trend, as it depends on three parameters (the asymptotic value b, the activation constant K act , and the Hill exponent p) which provide the flexibility required to tune up the sigmoidal shape of the representation.
The constitutive laws (20a) and (20c) are a mathematical translation of the adaptive response of the uveoscleral outflow pathway to an abnormal increase of the TM resistance which is a pathological precondition for the occurrence of glaucoma. 12Relation (20a) is a nonlinear Ohm's law representing the hydraulic component of the UncAR, while relation (20c) is a nonlinear drug-controlled current source representing the mechanochemical component of the UncAR. 25In the case of (20a), we see that if w(x) is negative, then L hyd Unc ¼ L Unc;b , in agreement with the physical consideration that no hydraulic AR is expected if the pressure is below its baseline value.In the case of (20c), we see that arbitrary increasing drug concentration does not significantly improve the efficacy of the mechanochemical AR (saturation effect).

VII. MODEL CALIBRATION AND PARAMETERS
In this section, we define the values attained by model parameters in baseline conditions, corresponding to a healthy physiological state of the eye of an individual.The baseline value of the AH volumetric flow rate of the drugdependent component of the unconventional pathway, Q drug Unc;b , is equal to zero because we assume that in healthy conditions there is no need to receive a drug treatment.

VIII. SOLUTION ALGORITHM
Replacing the constitutive laws (13) into Eq.(12d) and substituting the resulting expression for Q unc into Eq.( 16) yields where T m is the iteration function defined as According to Definition VIII.1, we see that the problem of finding the IOP x that solves the nonlinear equation ( 16) is equivalent to the problem of finding a fixed point of the function T m .To this purpose, we consider in the remainder of the article the following two situations: (i) increase of R Conv from baseline to pathological conditions; (ii) administration of a topical drug therapy to lower the increase of IOP resulting from (i).
The two above situations constitute the modular structure of the Virtual Computational Laboratory (VCL) that we have developed using the Matlab scientific environment. 27The mathematical definition and analysis of (i) and (ii) are addressed in Appendixes A and B, respectively, while their computational performance is investigated in Sec.IX.

IX. MODEL VALIDATION
In this section, we use the Virtual Computational Laboratory (VCL) proposed in this article to conduct a series of simulations according to the following protocol: In Phase 1, the VCL simulates the efficacy of the sole pressuredependent component of the UncAR to limit the increase of IOP due to a progressive increase of TM resistance, In Phase 2, the VCL simulates the efficacy of the drug-dependent component of the UncAR (in addition to the pressure-dependent component) to lower IOP by the progressive administration of a topical drug (Latanoprost 0.005%) according to a therapy of 14 days with 1 drop/day.
It is worth noting that the simulation of the drug administration therapy does not actually correspond to introducing in the model a "time" variable; rather, time is a parametric variation of the mass of administered drug according to 14 equal increments of 1 drop/day.In the model, the process of drug diffusion across the cornea and throughout the AC is not accounted for; instead, we assume that the quantity of drug increment that is actually adsorbed by the eye (i.e., the actual amount that reaches AC) is a percentage of 5% of the administered mass of drug increment.This choice is consistent with the design principles of ocular drug delivery systems. 29n the two-phase protocol, we set where j drug and j hyd are non-negative constants.In the remainder of the text, we set j drug ¼ 0:33, in such a way that (i) b drug unc satisfies the sufficient condition (B2f) and (ii) the predicted decrease of IOP falls within a range consistent with existing clinical data.Two choices are considered for j hyd , the first is j hyd ¼ 0, corresponding to a linear resistor model of the pressure-dependent component of the unconventional adaptive route, and the second is j hyd ¼ 0:99, corresponding to a nonlinear resistor model of the pressure-dependent component of the unconventional adaptive route.With both choices, the resulting value of b hyd unc satisfies the sufficient condition (B3b).Figures 11-17 illustrate the results of the two-phase simulation protocol.For each figure, the solid lines in magenta color are obtained with the linear resistor model, the solid lines in blue color are obtained with the nonlinear resistor model, and the dashed lines in black color represent the baseline value of the plotted quantity.The x-axis in the left panel represents the ratio g :¼ R Conv =R Conv;b between the TM resistance and its value at baseline conditions.Thus, a value of 3 on this axis means that the TM resistance is three times the value of the TM resistance at baseline.The x-axis in the right panel represents the number of days (14, in total) of the progressive administration of 1 drop/day of Latanoprost 0.005%.
The convergence test of the computational algorithm consisted of terminating the fixed-point iteration (B1a) at the first k Ã !0 such that jx ðk Ã þ1Þ À x ðk Ã Þ j x ðk Ã þ1Þ 10 À11 : For each simulated condition, k Ã did not exceed a value of 26.The total elapsed time was on the order of a few seconds, including figure plots, on a laptop equipped with Intel Core i5 processor.

F. Comments to AH mass conservation
In exact arithmetics, the computed value of Q AH should be exactly equal to Q b so that possible deviations from this value can be regarded as a measure of the accuracy of the model.The maximum absolute difference between Q b ¼ 2:75lLmin À1 and Q AH (computed with the approximate variable x obtained by using the nonlinear model of the unconventional hydraulic conductance) is 4:37 Â 10 À11 lLmin À1 , while the same absolute difference is equal to 2:2 Â 10 À15 lLmin À1 when using the linear model of the unconventional hydraulic conductance.These results demonstrate the high accuracy of the implemented algorithm in enforcing at the numerical level the property of water fluid mass conservation.

X. DISCUSSION
In this section, we analyze the physiological consistency of the simulation results illustrated in Sec.IX and the clinical relevance of the adaptive routes of the unconventional outflow pathway on the lowering of IOP.

A. The model for L hyd
Unc and the lowering of the IOP The left panel of Fig. 11 suggests that the IOP nonlinearly depends on the TM resistance R Conv even in the case where a linear resistor electric analog is used for the pressure-dependent component of the unconventional outflow pathway (corresponding to j hyd ¼ 0).To understand the cause of this behavior, we study the sensitivity S x ¼ S x ðgÞ of IOP (x) with respect to a relative change g in the TM resistance.We have where / b denotes the baseline value of /ðxÞ.The plot of S x ðgÞ in Fig. 16 reveals that (1) the IOP continuously increases since S x ðgÞ > 0; (2) the IOP increase becomes progressively smaller since S x ðgÞ is monotonically decreasing; (3) saturation of the IOP occurs only for unphysiologically elevated values of R Conv .
The above considerations support the conjecture that a linear model for the pressure-dependent adaptive route of the unconventional outflow pathway is not adequate to account for a significant lowering of the IOP as an adaptive response to a pathological increase of R Conv .To overcome this inadequacy, we, therefore, assume that j hyd > 0, in such a way that the hydraulic conductance L hyd Unc;b depends (nonlinearly) on the pressure drop p AC À p S , consistent with the conjecture of Costagliola et al. 13 The physiological motivation of the nonlinear model (20b) is based on the conjecture formulated by Liton and Gonzalez 11 of an adaptive regulatory mechanism in the cells of the outflow pathway capable to respond to the mechanical stress originated by changes in the IOP by triggering signals aimed at modulating the flow of AH.The sigmoidal shape of the Hill function ( 17) is an attempt to include a saturation of the IOP for elevated values of the TM resistance, in agreement with the concept of intraocular pressure homeostasis of keeping ocular pressure within relatively narrow acceptable bounds elaborated by Acott et al. 31 To test the efficacy of the nonlinear model (20b), we illustrate in Fig. 17 the model prediction of IOP as a function of g in the case where j hyd ¼ 10, corresponding to a value of b hyd unc that does not satisfy the sufficient condition (B3b).The solid line (blue color) in the left panel, obtained with the nonlinear model (20b), predicts a tendency to saturation of IOP to a value of about 19.2 mm Hg, to be compared with the maximum value of about 23 mm Hg predicted by the linear model.Since the usual criterion for the onset of glaucoma is an IOP higher than 21 mm Hg, 32 simulation predictions suggest that the nonlinear model for L hyd Unc is consistent with an adaptive response of the unconventional outflow pathway capable of restoring IOP homeostasis in the anterior segment of the eye. 31 The model for Q drug Unc and the lowering of the IOP The effect of drug administration on the uveoscleral outflow is subject to intensive investigation.Pilocarpine causes contraction of the ciliary muscle fibers and compression of extracellular space, thus reducing the uveoscleral outflow.
5][36][37] The experimental measurement of outflow through the unconventional pathway is quite a challenging task because of the difficulty in measuring the flow rate of AH as pointed out by Johnson et al. 15 (Sec.IV).This difficulty provides a motivation to devise alternate methods (indirect techniques) to gain information on the unconventional pathway outflow.In this perspective, the expression (20c) proposed in the present work can be regarded as an original instance of an indirect technique to model the AH volumetric flow rate across the UncAR, secondary to IOP-lowering drug administration.
In the simulations, we set F ¼ 0.05 (corresponding to an absorbed mass equal to 5% of the administered mass) and m onedrop ¼ 1:5625lg.Model predictions in the right panels of Figs.11 and 17 indicate that the reduction of IOP from the value corresponding to pathological conditions is equal to 7.19 mm Hg in the case of the linear resistor model, while the reduction in the case of the nonlinear resistor model is equal to 6.56 mm Hg if j hyd ¼ 0:99 and 3.48 mm Hg if j hyd ¼ 10.The predicted IOP reduction in the case j hyd ¼ 10 is compatible with the data reported by Wang et al. 38 where a mean reduction of 6:56 mm Hg was measured in glaucomatous patients treated with Latanoprost 0.005% for 4 weeks.

XI. CONCLUSIONS
The adoption of mathematical modeling and computational techniques to advance knowledge in ophthalmology is relatively recent but continuously increasing, as witnessed by the foundation of peerreviewed journals specifically devoted to the topic (cf.e.g., the website 39 ) and the development of editorial initiatives targeted to provide a unified framework to anatomy, physiology, and imaging. 40n this article, we propose and analyze a novel mechanism-driven mathematical model of aqueous humor dynamics in the human eye.The formulation is a compartment representation of the anterior segment of the eye based on the analogy between electric and fluid variables illustrated in Sacco et al. 9 The scheme is implemented in a Virtual Computational Laboratory that is used to investigate the adaptive routes exerted by the two components of AH drainage, the conventional and unconventional outflow pathways, to lower an elevated value of intraocular pressure in the eye due to an abnormal increase of the resistance of the TM as proposed by Costagliola et al. 13 The novelty of our model is the mathematical characterization of the unconventional outflow pathway by means of phenomenological constitutive laws that attempt to take into account the principal conjectures raised about the adaptive mechanisms that are put in action by the ciliary muscles to contrast an abnormal increase of the IOP.The proposed model of the UncAR is demonstrated: • to provide a correct physiological picture of the hydrodynamical conditions in the anterior segment of the eye; • to reproduce clinical observations in the cure of the eyes of patients affected by glaucoma.
The computational performance of the formulation is corroborated by an a priori theoretical analysis of the well-posedness of the fixed-point iteration which is used to numerically solve the nonlinear model equations.The importance of this analysis is that it provides quantitative information about the values of specific model parameters that are seen to play a relevant role in the function of the simulated system.
The proposed mathematical framework is affected by several limitations; specifically: (1) the electro-osmotic-oncotic pressure gradient which acts at the cellular and subcellular level to drive the fluid from the ciliary capillaries of the ciliary body into the PC is not accounted for in the description of the process of AH secretion; (2) the outflow pathway routes are described by means of phenomenological laws so that adaptive mechanisms occurring at the cellular scale level in both TM and ciliary muscle are not accounted for; (3) the processes of drug absorption across the cornea and drug diffusion into the AC of the eye are not included.
Regardless of the above limitations, the model proposed in this article has been demonstrated to be a useful tool to investigate the processes of circulation and drainage of aqueous humor in the anterior segment of the human eye.The obtained predictions of the intraocular pressure and AH flow are in agreement with ocular physiology and clinical data and provide an encouraging basis toward combining the proposed model with the modern techniques of Machine Learning and Artificial Intelligence proposed by Guidoboni et al. 41 and that promise to represent the winning approach to a routine implementation of "personalized medicine" in the cure of ocular diseases. 42

XII. FIGURES WITH SIMULATION RESULTS
In this section, we collect the figures illustrating the simulation results described in Sec.IX and discussed in Sec.X.

FIG. 1 .
FIG. 1. Inflow and outflow pathways of AH.Black arrows are used to indicate the secretion of AH from the ciliary body, its circulation from the PC, which is the region behind the iris and in front of the ciliary body and the lens, into the AC, and its drainage through the trabecular and uveoscleral outflow pathways.[Reprinted with permission from Ramakrishnan et al., Diagnosis and Management of Glaucoma, Aqueous Humor Dynamics (Jaypee Brothers Medical Publishers 2013), Chap.09.Copyright 2013 Authors, licensed under a Creative Commons Attribution (CC BY) license.]

FIG. 7 .
FIG. 7. Compartment equivalent representation of the anterior segment of the eye.
FIG. 8. Electric equivalent circuit representing AH circulation between the PC and AC and AH drainage through conventional and unconventional adaptive routes.The current source at the top of the scheme represents the AH volumetric flow rate that is secreted by the ciliary processes.The quantities L PC and L AC represent the hydraulic facilities of the PC and AC, respectively.The green arrow, representing the secreted AH volumetric flow rate, divides into two contributions, Q Conv , representing the portion of AH volumetric flow rate that is drained out through the conventional outflow pathway, and Q Unc , representing the portion of AH volumetric flow rate that is drained out through the unconventional outflow pathway.
12 and Q 23 are the AH volumetric flow rates between nodes 1 and 2 and nodes 2 and 3.The quantities Q Conv and Q Unc are the AH volumetric flow rates throughout the conventional and unconventional outflow pathways.The quantity Q hyd Unc is the hydraulicdependent component of the AH volumetric flow rate throughout the unconventional outflow pathway.The quantity Q drug Unc is the drugdependent component of the AH volumetric flow rate throughout the unconventional outflow pathway.
PC ; x IOP; y p AC and w ðy À p S Þ.Each equation in system(13) is the electric counterpart of the discrete analog of Darcy's law relating velocity to pressure gradient under the assumption that pressure is a linearly varying function of position.In particular, relations (13a)-(13c) are linear Ohm's laws, whereas relation (13d) is a nonlinear Ohm's law in which the hydraulic conductance L hyd Unc nonlinearly depends on the pressure drop w ¼ y À p S .
FIG. 9. Graph of the Hill function for b ¼ 1, in correspondence of several values of K act .

FIG. 10 .
FIG. 10.Blue line and blue y-axis: plot of Z as a function of the intraocular pressure x in the interval ½10; 30mmHg.Orange line and orange y-axis: plot of L hyd Unc as a function of Z.The mathematical expressions of Z and L hyd Unc can be found in Eqs.(19) and (20b), respectively.
et al.,26 the following definition holds.Definition VIII.1.Let G ¼ GðxÞ be a given continuous function.The real number f is a fixed point of G if f ¼ GðfÞ:

Phase 1 :
we solve (22a) assuming that (a) Q drug Unc ¼ 0; (b) the TM resistance is progressively increased from its baseline value R TM;b to a pathological value R TM;path ¼ 3R TM;b .The chosen value of R TM;path is consistent with the observation made in Ref. 28 and with the data of Ref. 18 which show a range between 0:31 6 0:18lL min À1 mmHg À1 (in human anterior segments treated with 100nM endothelin-1) and 0:49 6 0:26lL min À1 mmHg À1 (in control human anterior segments.)Phase 2: we solve (22a) assuming that: (a) the TM resistance is kept fixed to R TM;path ; (b) the mass of an IOP-lowering drug is progressively increased from its baseline value m b ¼ 0lg to a value m 14 ¼ 1:094lg, corresponding to a therapy of 14 days with 1 drop/day.The two-phase protocol has the following clinical interpretation.

FIG. 11 .
FIG. 11.Left panel: IOP as a function of g.Right panel: IOP as a function of the number of days of therapy.The red line is characterized by a minimum IOP of 15:73 mmHg, corresponding to a percentage variation of 31:37%.The green line is characterized by a minimum IOP of 15:36 mmHg, corresponding to a percentage variation of 33% (minimum drug efficiency reported in Ref. 30).The cyan line is characterized by a minimum IOP of 14:45 mmHg, corresponding to a percentage variation of 37% (maximum drug efficiency reported in Ref. 30).The values of the parameters for the drug-and pressure-dependent components of the UncAR are j drug ¼ 0:33 and j hyd ¼ 0:99.

FIG. 14 .
FIG. 14. Left panel: Q hyd Unc as a function of g.Right panel: Q hyd Unc as a function of the number of days of therapy.The values of the parameters for the drug-and pressure-dependent components of the UncAR are j drug ¼ 0:33 and j hyd ¼ 0:99.
The left panel corresponds to setting Q hyd drug ¼ 0. The right panel illustrates the drug-dependent component of the unconventional AH VFR as a function of the number of days of therapy (setting R Conv ¼ R Conv;path ): this quantity is not an outcome of model computation but a result of post-processing since Q drug Unc does not depend on the intraocular pressure x [see Eq. (20c)].

FIG. 17 .
FIG. 17. Left panel: IOP as a function of g.Right panel: IOP as a function of the number of days of therapy.The values of the parameters for the drug-and pressure-dependent components of the UncAR are j drug ¼ 0:33 and j hyd ¼ 10.
IOP path ¼ 21 mm Hg being the pathological value of the intraocular pressure, and we evaluate b hyd